/*
 *
 *                This source code is part of
 *
 *                 G   R   O   M   A   C   S
 *
 *          GROningen MAchine for Chemical Simulations
 *
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2009, The GROMACS development team,
 * check out http://www.gromacs.org for more information.

 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * If you want to redistribute modifications, please consider that
 * scientific software is very special. Version control is crucial -
 * bugs must be traceable. We will be happy to consider code for
 * inclusion in the official distribution, but derived work must not
 * be called official GROMACS. Details are found in the README & COPYING
 * files - if they are missing, get the official version at www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the papers on the package - you can find them in the top README file.
 *
 * For more info, check our website at http://www.gromacs.org
 */
#include "gmx_meancharges.h"



/*! \brief
 * Function that does the analysis for a single frame.
 *
 * It is called once for each frame.
 */
static int analyze_frame(t_topology *top, t_trxframe *fr, t_pbc *pbc,
              int nr, gmx_ana_selection_t *sel[], void *data)
{
    t_analysisdata     *d = (t_analysisdata *)data;

    uint nx=(uint) ngrid[0];
    uint ny=(uint) ngrid[1];
    uint nz=(uint) ngrid[2];
/*
    Charges=malloc(sizeof(real)*top.atoms.nr);

	for(i=0;i<top.atoms.nr;i++) {
		Charges[i]=top.atoms.atom[i].q;
	}
*/
	if (bCenter) center_x(ecenter,fr->x,fr->box,fr->natoms,ncent,cindex);
	t_pbc pbc;

	set_pbc(&pbc,ePBC,fr->box);
	if (bRmPBC)
		rm_pbc(&(top->idef),ePBC,fr->natoms,fr->box,fr->x,fr->x);

	if (bPBCcomAtom) {
		switch (unitcell_enum) {
			case euRect:
			put_atoms_in_box(fr.box,fr.natoms,fr.x);
			break;
			case euTric:
			put_atoms_in_triclinic_unitcell(ecenter,fr.box,fr.natoms,fr.x);
			break;
			case euCompact:
			warn = put_atoms_in_compact_unitcell(ePBC,ecenter,fr.box,
					fr.natoms,fr.x);
			break;
		}
	}
	if (bPBCcomRes) {
		put_residue_com_in_box(unitcell_enum,ecenter,
				fr.natoms,atoms->atom,ePBC,fr.box,fr.x);
	}
	if (bPBCcomMol) {
		put_molecule_com_in_box(unitcell_enum,ecenter,
				&top.mols,
				fr.natoms,atoms->atom,ePBC,fr.box,fr.x);
	}

	if(bTrans) {
		for(int i=0; i<fr.natoms; i++)
			rvec_inc(fr.x[i],trans);
	}

	if (bOXT) {
		acc_x(ncent, cindex, fr.x, x);
	}

    /* Here, you can do whatever analysis your program requires for a frame. */
    printf("%10.3f", fr->time);

    printf("\n");

    /* We need to return 0 to tell that everything went OK */
    return 0;
}

/*! \brief
 * Function that implements the analysis tool.
 *
 * Following the style of Gromacs analysis tools, this function is called
 * \p gmx_something.
 */
int main(int argc, char *argv[])
{
    const char         *desc[] = {
        "This is a template for writing your own analysis tools for",
        "Gromacs. The advantage of using Gromacs for this is that you",
        "have access to all information in the topology, and your",
        "program will be able to handle all types of coordinates and",
        "trajectory files supported by Gromacs. In addition,",
        "you get a lot of functionality for free from the trajectory",
        "analysis library, including support for flexible dynamic",
        "selections. Go ahead an try it![PAR]",
        "To get started with implementing your own analysis program,",
        "follow the instructions in the README file provided.",
        "This template implements a simple analysis programs that calculates",
        "average distances from the a reference group to one or more",
        "analysis groups.",
    };

		/* Command-line arguments */


    t_pargs             pa[] = {
    		{"-grid", FALSE, etRVEC, {&ngrid}, "Grid size"},
            { "-pbc", FALSE, etENUM,
            		{ pbc_opt },
                "PBC treatment (see help text for full description)" },
            { "-ur", FALSE, etENUM,
                { unitcell_opt }, "Unit-cell representation" },
            { "-dipole", FALSE, etBOOL,
                { &bDipole }, "The center molecule is a dipole" },
            { "-filter", FALSE, etREAL,
                {&filter}, "The center molecule is a dipole" },
            { "-xplor", FALSE, etBOOL,
                {&bXplor}, "Write a potential map in Xplor format" },
            { "-trans", FALSE, etRVEC,
                  { trans },
                  "All coordinates will be translated by trans. This "
                  "can advantageously be combined with -pbc mol -ur "
                  "compact." },
            { "-boxcenter", FALSE, etENUM,
                { center_opt }, "Center for -pbc and -center" },
    };

    /* The second argument is for demonstration purposes only */

      t_filenm fnm[] = {
    { efTPS,  NULL,  NULL, ffREAD },   /* this is for the topology */
    { efTRX, "-f", NULL, ffREAD },      /* and this for the trajectory */
    { efXVG, "-dist", "distr",    ffOPTWR  },
    { efTRX, "-oxt","coord.xtc", ffOPTWR },
    { efNDX, NULL,  NULL,         ffOPTRD },
  };
#define NFILE asize(fnm)
#define NPA asize(pa)

    /* Output files */


    gmx_ana_traj_t       *trj;
    output_env_t          oenv;
    t_analysisdata        d;
    int                   ngrps;
    gmx_ana_selection_t **sel;

    CopyRight(stderr, argv[0]);
    /* Here, we can use flags to specify requirements for the selections and/or
     * other features of the library. */
    gmx_ana_traj_create(&trj, ANA_REQUIRE_TOP);
    gmx_ana_set_nrefgrps(trj, 1);
    gmx_ana_set_nanagrps(trj, 1);
    /* If required, other functions can also be used to configure the library
     * before calling parse_trjana_args(). */
    parse_trjana_args(trj, &argc, argv, PCA_CAN_VIEW,
                      NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
                      &oenv);


	pbc_enum   = nenum(pbc_opt);
    bPBCWhole  = pbc_enum==epWhole;
    bPBCcomRes = pbc_enum==epComRes;
    bPBCcomMol = pbc_enum==epComMol;
    bPBCcomAtom= pbc_enum==epComAtom;
    bNoJump    = pbc_enum==epNojump;
    bCluster   = pbc_enum==epCluster;
    bPBC       = pbc_enum!=epNone;
    bTrans     = opt2parg_bSet("-trans",NPA,pa);
    bFilter    = opt2parg_bSet("-filter",NPA,pa);
    unitcell_enum = nenum(unitcell_opt);
    ecenter    = nenum(center_opt) - ecTric;

    bRmPBC = bPBCWhole || bPBCcomRes || bPBCcomMol;

    distfile = opt2fn_null("-dist",NFILE,fnm);
    bOXT = opt2bSet("-oxt",NFILE,fnm);

    if(distfile) fp = xvgropen(distfile,"","","",oenv);
    t_trxstatus *status_out=NULL;
	if (bOXT) {
		flagsb = flags | TRX_READ_X;
		status_out = open_trx(opt2fn("-oxt",NFILE,fnm),"w");
	}

	bIndex = ftp2bSet(efNDX, NFILE, fnm);

//	if (bIndex)
//		printf(" Success !!\n");

	rvec       *xtop;

	read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, NULL,
			box, TRUE);
	sfree(xtop);
	t_atoms      *atoms=NULL;
	atoms = &top.atoms;

	if(bOXT){
		int natom=atoms->nr;
		x=calloc(natom,sizeof(rvec));

	}
	if (bCenter) {
		if (bIndex) {
			printf("Select group for centering\n");
			rd_index(ftp2fn_null(efNDX, NFILE, fnm), 1, &ncent, &cindex, &grpnm);
		} else {
			gmx_fatal(FARGS, "Cannot do centering without Index file. Abort.");
		}
	}


    /* You can now do any initialization you wish, using the information
     * from the trj structure.
     * In particular, you should store any command-line parameter values that
     * the analysis part requires into d for them to be accessible. */

    /* First, we get some selection information from the structure */
    gmx_ana_get_refsel(trj, 0, &d.refsel);
    gmx_ana_get_nanagrps(trj, &ngrps);
    gmx_ana_get_anagrps(trj, &sel);



    /* Now, we do the actual analysis */

    gmx_ana_do(trj, 0, &analyze_frame, &d);

    /* Now, the analysis has been done for all frames, and you can access the
     * results in d. Here, you should post-process your data and write out any
     * averaged properties. */



    /* Here, we could free some memory, but this usually not necessary as we
     * are going to quit anyways. */

    thanx(stderr);
    return 0;
}

